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Field of the Invention 
The field of the invention relates to blood 
circulation in the human body and more particularly, to 
blood flow within the human brain. 

10 Description of the Invention 

Stroke is a leading cause of death and disability 
in the United States, with significant socioeconomic 
impact. Therapeutic options for occlusive 
cerebrovascular diseases or trauma include a variety of 

15 reconstructive procedures such as endarterectomy , vessel 
transposition, bypass, angioplasty and thrombolysis; all 
of which share the common goal of enhancing the cerebral 
circulation . 

However, because of the many individual 

2 0 variabilities between patients, assessment of the 

potential merit of such procedures is difficult and has 
mainly relied on clinical trials involving large groups 
of patients over long periods of time and at staggering 
costs (e.g., the NIH EC-IC bypass and carotid 

2 5 endarterectomy studies) . The ethical dilemma of 

randomizing a patient to an ineffective or even harmful 
therapy not withstanding, selection of the optimal 
cerebral vasdular reconstructive procedure in a given 
patient with complex occlusive patterns remains 

3 0 difficult and has to date relied on an intuitive 

process . 
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The ability to reliably simulate such procedures 
for any particular patient, increases the likelihood of 
selecting the optimal procedure. This invention applies 
a computerized model of the cerebral circulation, to 
patients with cerebrovascular diseases to simulate any 
of a number of procedures, thereby providing a tool for 
selecting the optimal procedure for any given patient. 
The benefit of this invention provides improved patient 
outcome, significant cost savings and major impact on 
the way cerebrovascular reconstructive. procedures are 
performed . 

Charbel et al . , 2"*^ International Skull Base 
Congress , 7^^ Annual Meeting of the North American Skull 
Base Society, (San Diego, July, 1996) discussed the 
applications of an earlier computerized model of the 
cerebral circulation in skullbase surgery. The role of 
computerized modeling in cerebrovascular surgery was 
discussed by Charbel et al . at the 11^^ International 
Congress of Neurological Surgery (Amsterdam, July 1997). 
Charbel et al . discussed the use of a computerized 
predictor for the tolerance of carotid artery sacrifice 
at 11^^ International Congress of Neurological Surgery 
(Amsterdam, July 1997). Applications of modeling to 
cerebral revascularization are also contemplated. 

Summary 

A method and apparatus are provided for modeling 
circulation in a living subject. The method includes 
the steps of developing a model for living subjects in 
general and correcting the model to substantially 
conform to the overall cerebral physiology of the living 
subject. The method further includes the step of 



calculating a cerebral flow of the living subject based 
upon the corrected model and a selected cerebral blood 
flow intervention (such as a bypass or carotid 
sacrifice) . 

Brief Description of the Drawings 
FIG. 1 is a block diagram of a modeling system for 

cerebral circulation in accordance with an illustrated 

embodiment of the invention; 

FIG. 2 is a cerebral circulation portion of the 

model used by the system of FIG. 1; 

FIG. 3 is an operator interface used by the system 

of FIG. 1; 

FIG. 4 is a system calibration screen used by the 
operator interface of FIG. 3; 

FIG. 5 is a vessel selection screen used by the 
system of FIG. 1; 

FIG. 6 is a vessel size measurement screen used by 
the system of FIG. 1; 

FIG. 7 depicts a series of images generated by the 
system of FIG. 1; 

FIG. 8 depicts a histogram of an MIP image for 
determining an iso-surface by the system of FIG. 1; 

FIG. 9 depicts a sagittal plane and vasculature 
obtained using the Marching-Cube Algorithm by the system 
of FIG. 1; 

FIG. 10 depicts a screen used in flow measurement 
by the system of FIG. 1; 

FIG. 11 depicts color-coded flow of a region of 
interest used by the system of FIG. 1; 

FIG. 12 depicts a graphical representation of blood 
flow used by the system of FIG. 1; 



FIG. 13 depicts a vasculature structure and slice 
image used by the system of FIG. 1; 

FIG. 14 depicts a flow phantom and surface 
rendering of the flow phantom used by the system of FIG. 
5 1; 

FIG. 15 depicts slice sections of the flow phantom 
of FIG. 14; 

FIG. 16 depicts Cine PCMR magnitude images used by 
the system of FIG. 1; 
10 FIG. 17 depicts a comparison of flow measurements 

used by the system of FIG. 1; 

FIG. 18 depicts a comparison of flow measurements 
for a pulsatile flow phantom used by the system of FIG. 
1; 

15 FIG. 19 depicts 3D surface renderings of a 

particular cut and three misalignment angles used by the 
system of FIG. 1; 

FIG. 20 depicts a 3D localization of 4 particular 
intercranial arteries provided by the system of FIG. 1; 
20 FIG. 21 depicts a generic 108-vessel cardiovascular 

model used by the system of FIG. 1; 

FIG. 22 shows definition sketches of the model used 
by the system of FIG. 1; 

FIG. 23 depicts a seven-sector cerebral model used 
25 by the system of FIG. 1; 

FIG. 24 is a flow chart of the sector model used by 
the system of FIG. 1; 

FIG. 25 depicts an angiogram, 3D localizer image 
and flow waveform of a patient provided by the system of 
30 FIG. 1; 

FIG. 2 6 depicts case results produced by the system 
of FIG. 1; 
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FIG. 27 depicts further case results produced by 
the system of FIG. 1; 

FIG. 2 8 depicts simulated CBF distributions and BOT 
of Case 2 performed by the system of FIG. 1; and 

FIG. 29 is a flow chart for sequencing calculations 
for the model used by the system of FIG. 1. 

Detailed Description of an Illustrated Embodiment 
The circulation around the "Circle of Willis" as 
first described in 1664 was initially thought to be a 
poorly functioning anastomotic network at the base of 
the brain. Willis, Annals of Medical History 2:81-94 
(1940) . Later it became accepted that it acts as the 
main distribution center for cerebral blood flow. In 
health, it distributes blood proportionately to each 
part of the brain; in disease, when the blood supply is 
decreased or diminished, it can redistribute flow in an 
effort to maintain cerebral hemodynamic homeostasis. 

Because the cerebral circulation network is encased 
by the skull, it is very difficult to measure flow and 
pressure in its blood vessel constituents in a direct 
manner. Therefore computer simulation using the system 
10 becomes an attractive way to predict flows and 
pressures in the cerebral circulation network. Fluid 
dynamic-based computer simulation offers the convenience 
of predicting pressures and flows at almost any desired 
section in the circulating system. Moreover, it not 
only can be used to estimate the flow and pressure under 
health and disease situations, but also can be used to 
predict the result of treatment procedures. 

The simulation of cerebral circulation presents a 
range of challenging fluid dynamic problems, including: 



modeling the non-Newtonian properties of blood; dealing 
with "physiological" unsteady pulsatile flow; modeling 
the elasticity of vessel walls; and modeling moving 
boundaries caused by vessel wall elasticity. In 
addition, because the cerebral circulation network is an 
interconnected three-dimensional arterial network, the 
question arises of how the curvature of the artery is 
modeled. The asymmetric and three-dimensional 
characteristics of bifurcations are also very important 
issues. 

Computer simulation started with models of the 
dog's cerebral circulation system. Clark et al . , Acta 
Neurol. Scandinav . , 43:189-204 (1967), built a computer 
model for one-dimensional, linear, steady laminar flow 
and compared the result with an engineering model built 
by the same group [Himwich et al . , Archive of Neurology, 
13:164-172 (1965); Himwich et al . , Archive of Neurology, 
13:173-182 (1965)]. Cooper, Ph.D. thesis, Washington 
University (1970), studied the pulsatile flow effect. 
Chao & Hwang, TlT Journal of Life Science, 2:3-81 (1972) 
simulated non-linear pulsatile flow in their model. 
Himwich & Clark, Pathology of Cerebral Microcirculation , 
J. Cervos -Navarro , ed. (Walter de Gruyter, New York: 
1974), pp. 140-152, and Clark & Kufahl, Proc . of 1^^ 
Intl. Conf . Cardiovascular System Dynamics , MIT Press 
(Boston, 1978), pp. 380-390, simulated the pulsatile 
flow in distensible vessels. Hillen et al . , J". 
Biomechanics, 15:441-448 (1982), developed a non-linear, 
one-dimensional model, simplified to only five vessels, 
that accounted for pulsatile flow and elastic vessel . 
walls. Kufahl & Clark, J". Blomecha.nlca.1 Engineering, 
107:112-122 (1985), developed a one-dimensional finite- 



difference model with distensible vessels. Kufahl & 
Clark's dog model contained thirty-five vessels and 
allowed both steady flow and pulsatile flow to be 
studied . 

5 Beginning in the late eighties, investigators began 

to simulate the human cerebral circulation. The process 
of adapting the animal models to generally simulate 
human cerebral circulation was not trivial . 

The cerebral circulation network of a human is more 

10 complex than that of an animal. Considerations to be 
addressed included the number of arteries it was 
necessary to simulate, the selection of those arteries, 
and the configuration of the chosen arteries to 
represent the functionality of the Circle of Willis. 

15 Cerebral circulation in humans is not consistent 

among individuals. The number of arteries, and 
parameters such as artery length and diameter differ 
from person to person, limiting the utility of standard 
data for individual cases. A successful model must be 

2 0 very flexible and robust to handle various simulation 
conditions . 

The difficulty of obtaining direct measurements in 
humans limits the ability to accumulate parametric 
information and to refine the model . 

25 Clark et al . , Neurological Res., 11:217-230 (1989), 

developed an elaborate model for the human cerebral 
circulation network following the mathematical method 
used by Kufahl & Clark, J". Biomecha.nica.1 Engineering , 
107:112-122 (1985). The new model included seventy- 

30 three vessels representing the basic Circle of Willis. 
FIG. 2 is the schematic drawing of the 73 vessel model 
of Clark et al . After adding naturally occurring 




anastomoses, the total niimber of vessels increased to 
eighty-five. When artificial anastomoses were imposed, 
the number of vessels increased to eighty-seven. 

Duros et al . , Neurological Res,, 13:217-223 (1991), 
5 built a model that not only contained the cerebral 

arteries but also the human body main supply arteries. 
The Duros et al . model was used to simulate the rupture 
condition of aneurysm. 

Beside the fluid dynamic models, several electrical 
10 models were built based on the similarity of the 

governing equations of electrical circuits and one- 
dimensional linear flow. Moreover, electrical networks 
v3 are good at simulating networks with capacitance and 

^ resistance. Electrical network models are also well 

15 understood and provide convenient abstractions. Roller 
ffl Sc Clark, J, Biomechanics, 2:244-251 (1969) simulated the 

J' pulsatile flow and flexible vessel wall using 

O transmission line theory. Hillen, Comput , Biol. Med., 

^ 24:103-118 (1994), built an electrical model using 

•sew. 

20 transmission line equations similar to the linearized 
Navier-Stokes equation and vessel wall deformation 
equation . 

All existing cerebral circulation models are one- 
dimensional. All fluid dynamic-based models begin with 

25 the same set of governing equations, namely the one- 
dimensional continuity equation and Navier-Stokes 
equation. In those models, the following assumptions 
must be made in order to derive the equations [Long et 
al., J". Fluid Mech, , 55:493-511 (1972); Imaeda, Ph.D. 

30 thesis, Univ. of Waterloo (1975); Imaeda et al . , J". 
Biomechanics, 13:1007-1021 (1980); Nerem et al . , 
Handbook of Engineering (McGraw-Hill, New York: 1987), 
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pp. 21.1-21.21]: (i) the flow is radially symmetric and 
includes well -developed laminar flow; (ii) blood is an 
incompressible, Newtonian fluid; and (iii) external 
forces (e.g., gravitational force) are negligible. 

Because the most important components of interest 
are. flow and pressure distributions, the one-dimensional 
integrated continuity and momentum equations are more 
convenient to use . The integrated governing equations 
are [Kufahl & Clark, J. Biomechanical Enginoering , 
107:112-122 (1985); Raines et al . , Proc . of the Summer 
Computer Simulation Conference (MIT Press, Boston: 
1975), pp. 890-900; Raines et al . , J. Biomechanics, 
7:77-91 (1974); Proenta et al . , J. Biomechanical Eng., 
108:161-167 (1986) ] : 

continuity equation: 



(1) 



dx dt 



20 



momentuin elation: 



(2) 



dt dx 



A dP ^ TCDt 
P dx p 
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where Q is a quantum flow and A is cross- sectional area. 

t is the shear stress at the wall. The convective term 
contributes to the non-linearity of the governing 
equation set. 
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The main differences between the various models are 
how they define the problem (e.g., whether the flow is 
steady or pulsatile, whether the vessel wall is rigid or 
elastic, and whether a linear or non-linear governing 
equation is used to describe the flow field) and how 
they solve the problem (e.g. whether an analytical 
method or numerical method is used) . 

Even when the assumption is made that blood is 
Newtonian and flow is laminar, the governing equations 
are still non-linear because of the existence of the 
convective term in the momentum equation. Other non- 
linearities are included in the equation of state. For 
pipe flow, the convective term can be ignored because 
the gradient of the velocity is very small. The result 
is a linearized Navier-Stokes equation. If it is 
further assiimed that the vessel wall is rigid and the 
flow is steady, the problem becomes Poiseuille flow, 
which can be solved analytically by applying the well- 
known Hagen-Poiseuille formula. 

(3) Q^^rR!k_z_^ 

where R is the pipe radius and £ is the pipe length. 

For the rigid vessel unsteady pulsatile flow, 
Wormersley, Phys . Med. Biology, 2:178-187 (1957), 
provided an analytical solution for a single vessel. 
Ling & Attack, J, Fluid Mech, , 55:493-511 (1972), 
extended the Wormersley model by taking into account the 
consideration of non-linear effects and distensibility 
of the blood vessel. The finite-difference technique 
was employed to solve the problem numerically. To date, 
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most existing numerical models can handle the non- 
linear, pulsatile flow and distensible vessel wall 
problem. 

From the definition of the problem, computer 
modeling can be divided into two categories. 

The first type of models applied the linear 
governing equation and rigid vessel assumption. For 
steady flow, the solution is given by the Hagen- • 
Poiseuille formula. For unsteady flow, the solution is 
given by the Wormersley equations. Examples of models 
of this type include Clark et al . Acta Neurol. 
Scandinav., 43:189-204 (1967), and Hillen et al . , J. 
Biomechanics, 21:807-814 (1988). Because some of the 
electrical models [Roller and Clark, J. Biomechanics , 
2:244-251 (1969); Hillen, Comput . Bio, Med., 24:103-118 
(1994)] were derived using the Hagen-Poiseuille formula 
as a model, they belong to this type. The advantages of 
a model with linear governing equation and rigid vessel 
assumptions is that the model is simple and an 
analytical solution exists. Thus, sensitivity analysis 
of parameters can be performed. It also offers the 
possibility of gaining insight into mechanisms that 
govern the flow in the cerebral network. 

The other type is non-linear, pulsatile flow and 
distensible vessel models. Examples of models of this 
type include Kufahl & Clark [Kufahl and Clark, J. 
Biomechanical Engineering, 107:112-122 (1985); Clark et 
al.. Neurological i^esearch, 11 : 217-230 (1989); Kufahl, 
Ph.D. thesis (Univ. oflllinois, Urbana: 1980)] Hillen 
et al . [Hillen et al . , J. Biomechanics, 15:441-448 
(1982); Hillen et al . , J. Biomechanics, 19:187-194 




(1986)] and Duros et al . , Neurological Research, 13:217- 
223 (1991) , 

The advantage of this kind of model is that the 
non-linear, pulsatile flow and distensible vessels are 
5 more accurate than the first type in describing the 
physical behavior of blood flow in arteries. For 
instance. Ling and Attack, J. Fluid Mechanics , 55:493- 
411 (1972), point out that the non-linear term cannot be 
neglected. Kufahl and Clark, J. Biomechanical 
10 Engineering, 107:112-122 (1985), showed that pulsatile 

flow resulted in a 10 percent increase of the total flow 
^ and contributed to some flow redistribution. 

For the first type of problem, analytical solutions 
X exist. For the second type of problem, numerical 

P 15 techniques must be used. The most commonly used 

s _ = 

^ computational fluid dynamics tools are finite-difference 

Ul (FD) , finite-element (FE) and finite volume (FV) . 

p In the task of simulating the human cerebral 

5=^ circulation system, the finite-difference technique is 

p 2 0 widely used [Kufahl & Clark, J". Biomechanical 

^2 Engineering, 107:112-122 (1985); Clark et al . , 

Neurological Research, 111:217-230 (1989); Duros et al . , 
Neurological Research, 13:217-223 (1991); Kufahl, Ph.D. 
thesis (Univ. of Illinois, Urbana, 1980)] because the 
25 application of finite-difference approach is 

straightforward and it can achieve certain accuracy with 
modest computational resources. However the finite- 
difference approach requires high regularity of the grid 
used, and this restrains the approach to be used in 
3 0 solving complex geometry problems. 

In contrast, the finite-element approach can handle 
the complex geometry problem very easily, but the 
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shortcoming of this method is that it needs much more 
computation time than the finite-difference method. 

Cerebral network models usually contain hundreds of 
arteries and tens of bifurcations, which makes the use 
5 of the finite-element approach difficult using presently 
available computer equipment. To apply the finite 
element method in the network simulation, powerful 
computers are needed. With advances in computation 
technology, the finite-element method can be applied in 
10 cerebral network simulation. The present invention 

currently incorporates the finite-difference approach. 

Stroke is a leading cause of death in the United 
States, as well as in most western countries. Among the 
2r many approaches available for cerebral revascularization 

p 15 (e.g., angioplasty, endarterectomy , bypass and 

if embolectomy) , the procedure of choice for each 

yi particular patient is, at least in theory, the one that 

restores cerebral blood flow. However because of the 
ip complex architecture of cerebral circulation, each 

2 0 patient has a unique dimensional structure. As a 
result, any surgical treatment can have different 
effects on different patients. 

Computer models can simulate the cerebral 
circulation under "normal" conditions. Computer models 
25 can also be used to predict the results of potential 

treatment procedures. However, there is no other system 
currently available which forms a comprehensive model 
customized to the patient or which allows a user to 
perturb that model at will . 
30 Hillen et al., J. Biomechanics, 15:441-448 (1982), 

built a non-linear one-dimensional model to study the 
functional significance of the Circle of Willis. The 




model consisted of two afferent and two efferent 
arteries connected by the posterior communicating 
artery. Hillen, et al . found that in normal cases, the 
flow in the posterior communicating artery was towards 
5 the posterior cerebral artery, and that the flow 

direction in posterior communicating artery depended on 
the ratio of peripheral resistance. 

Kufahl and Clark, J. Biomechanical Engineering, 
107:112-122 (1985), developed a 35-vessel computer model 
10 for dog circulation. The model was non-linear, one- 
dimensional with pulsatile flow and distensible vessel 
^ walls. The computer simulation revealed the large drop 

^ in pressure over the length of the middle cerebral due 

j to the large fluid friction at the vessel wall. Large 

□ 15 pulses of flow were found in the common carotid but 

^ smaller pulses occurred farther downstream, 

yi The comparison by Kufahl and Clark of running the 

model with steady flow versus pulsatile flow showed that 
^ the total flow rate increased because of the pulsatile 

2 0 flow. Kufahl and Clark also found flow redistribution 
C in some arteries . 

A series of stenoses were simulated by Kufahl and 
Clark to study the efficacy of the Circle of Willis in 
maintaining the flow under arterial disease. The middle 
25 cerebral flow was found to be well-maintained under 
different cases. 

Hillen et al . , J. Biomechanics, 21:807-814 (1988), 
extended their previous model [Hillen et al . , J". 
Biomechanics 15:441-448 (1982)] to include 18 vessels. 

3 0 The effect of asymmetry on the flow distribution in 

afferent and efferent vessels and the possibility of the 
occurrence of a "dead point" in the posterior 
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communicating artery were studied. The authors found 
that a slight asymmetric change in the right posterior 
communicating artery (doubling the diameter) could 
result in noticeable effects in the flow distribution 
(flow in the left carotid artery exceeded that in the 
right carotid artery and a small flow was present in the 
anterior communicating artery) , but the pressure seemed 
to be unaffected. The authors pointed out that a dead 
point in the posterior communicating artery cannot occur 
in humans . 

Hillen et al . , J. Biomechanics, 21:807-814 (1988), 
simplified their model by ignoring pulsatility and 
vessel wall elasticity. From the comparison of the 
simple model and the previous model, they concluded that 
the deletion of pulsatility and vessel wall elasticity 
caused only a slight decrease of the flows without any 
flow redistribution. This finding was in contrast to 
the study of Kufahl and Clark, J". Biomechanical 
Engineering, 107:112-122 (1985), in which the flow 
redistribution was found to be due to pulsatile flow. 
The authors concluded the contradiction was caused by 
the ways in which the frictional resistances were 
determined in the different models. Hillen et al . used 
the Hagen-Poiseuille resistance, whereas Kufahl and 
Clark calculated wall friction from a quasi-two- 
dimensional velocity profile that is approximated by a 
sixth degree polynomial. 

Clark et al . , Neurological Research, 11:217-230 
(1989), built four computer models for the human Circle 
of Willis. The first model contained 73 vessels ("model 
73") representing the basic circle and served as a 
benchmark. The second model was constructed by adding 
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the naturally-occurring secondary anastomoses, the total 
number of vessels increased to 85 ("model 85"). Adding 
an artificial anastomosis between the frontal artery and 
the middle cerebral artery in model 73 and model 85 
5 resulted in model 7 5 and model 88, respectively. 

In order to evaluate the ability of each of the 
anastomotic vessels under normal conditions, Clark et 
al . studied five cases. In each case, a slight increase 
in the diameter of the anastomoses was observed over the 
10 previous model. The results showed that the EC-IC 

bypass was very beneficial when natural anastomoses were 
absent . 

I z 

^ In another five cases, a 90-percent stenosis was 

^ imposed in the middle cerebral artery. The same studies 

Q 15 were performed by Clark et al , as before. They found 

that the anastomotic channels helped supply the ischemic 
[fl areas, the larger the diameter, the better the supply. 

1^ Again, the artificial anastomoses were most beneficial 

,B when the natural anastomoses were not present. 

Iz 20 Later, Charbel et al . ["Predictive value of a 

computerized model of the cerebral circulation" , 44^^ 
Annual Meeting- of Congress of Neurological Surgeons 
(Chicago, 1994); "Validation and clinical potential of a 
computerized model of the cerebral circulation" , 1^^ 
25 Annual Meeting on the Joint Section on Cerebrovascular 
Surgery of the AANS & CMS (1996)] introduced both 
quantitative, semi-quantitative, as well as direct 
quantitative methods of validating a patient-specific 
computer model by directly measuring flow in blood 
3 0 vessels during surgery and comparing the values with the 
computer model . More recently, the same group began 
utilizing magnetic resonance as another valuable 
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validation tool, "Phase Contrast MR flow measurement 
system using volumetric flow constrained image 
interpolation and color coded image visualization", 47th 
Annual Meeting of Congress of Neurological Surgeons, 
5 (New Orleans, 1997), thus further bringing computer 
modeling within the reach of the clinician. 

Duros et al . , Neurological Research, 13:217-223 
(1991) , built a 69-vessel model to simulate the cerebral 
circulation with an aneurysm. The mathematical method 
10 used in the model was the same as Kufahl and Clark, J". 
Biomechanical Engineering, 107:112-122 (1985). The 

^ Duros et al . model differed from the other models in 

U 

^ that the model not only contained cerebral arteries, but 

^ also contained 3 0 main supply arteries to different 

O 15 organs of the human body. 

p The aneurysm was balloon- shaped, elastically 

U1 tapered with zero distal flow and was placed at the 

Q junction of the internal carotid, the anterior cerebral 

tfJ and the middle cerebral arteries. To simulate the 

2 0 condition where a rupture may happen, several parameters 
C were adjusted: all terminal vessels' resistance values 

^ were enlarged by a factor of 10 and two 80 percent 

stenoses were placed in the middle cerebral artery and 
anterior cerebral artery to flank the aneurysm. 
25 Systemic pressure was increased to 150 mm Hg to 

represent hypertension. The compliance coefficient of 
the aneurysm was set to 1.5 to represent a stiff wall 
condition. 

Duros et al . focused on the pressure propagation 

3 0 inside the aneurysm and found that the pressure did not 

change with the neck diameter, but the pressure peak 
value increased by increasing the sack diameter. In 
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order to achieve a high pressure (310 mm Hg) indicative 
of hypertension which may trigger the rupture of the 
aneurysm, an increased number of reflecting sites in 
both the near and far fields, and arteriosclerotic 
5 arteries were needed. 

Flow branching at arterial bifurcations is an 
important factor in the progression of vascular disease. 
This blood flow through the artery bifurcation was 
extensively studied both by experimental and numerical 
10 methods. 

Liepsch, Biorheology, 21:571-586 (1984), studied 
the non-Newtonian fluid in a T-shaped bifurcation. 
5 Perktold et al . , Biorheology, 26:1011-1030 (1989), 

]£ Studied a Y-shaped bifurcation with an aneurysm in the 

O 15 cerebrovascular tree. Walburn, J". Biomechanical 

55 Engineering, 104:66-88 (1982), studied steady flow and 

pulsatile flow through the aortic bifurcation, secondary 
Q flow was not observed during the pulsatile flow but was 

-S observed during the steady flow. Fukushima et al . , J. 

H 20 Biomechanica.1 Engineering, 110:161-171 (1988), found 

^ secondary flow in both steady and pulsatile flow, 

however, the influence of the secondary flow is smaller 
during the pulsatile flow than during the steady flow. 
To take the importance of the secondary flow into 
25 consideration, a few three-dimensional computer models 
have been built. Wille, J. Biomechanical Engineering, 
6:49-55 (1984), built a single three-dimensional model 
of aortic bifurcation with steady flow. Yung et al . , J". 
Biomechanical Engineering, 112:189-297 (1990), studied 
3 0 the steady flow in the aortic bifurcation. In order to 
generate the mesh conveniently, Yung et al . used a non- 
physiological area ratio of 2.0 to produce flow 
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geometry. Perktold & Peter, J. Bioinecha.nica.1 
Engineering, 12:2-12 (1990) studied the wall shear 
stress in a three-dimensional model of a T-shaped 
bifurcation. Perktold [J". Biomechanical Engineering, 
5 13:464-475 (1991); J". Biomechanics, 24:409-420 (1991)] 
analyzed non-Newtonian characteristics, wall shear 
stress and pulsatile flow in a three-dimensional model 
of carotid bifurcation. Rindt et al . , built a three- 
dimensional model for carotid bifurcation, but they used 

10 steady flow and rigid vessel wall. 

The simplest blood vessel wall model is a rigid 
tube. However, vessel wall elasticity has an important 
effect on the blood flow wave propagation. For a large 
artery, the vessel wall is composed of three distinct 

15 layers, an intima, a media and an adventitia. Each of 
these layers has a unique function. The vessel wall is 
not only elastic but also viscoelastic . The vessel 
elasticity was modeled by a relationship between blood 
vessel cross-sectional area and local pressure. Raines 

2 0 proposed a widely used model: 

( 4 ) A{p,x)= A{p„ ,x)+pln{p/pj 

Another quadratic form used by Porenta, Balar, and 
25 Stergiopulos is : 

( 5 ) A{x) = A„ {xi + cAp-Po)+CAp- Po T J 

where (3, C'o and C'l are constants which are determined 

3 0 by experiment. 
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These models are purely elastic models. However 
Patel and Vaishnav verified the existence of the 
viscoelasticity in arterial walls through a dynamic 
experiment. Reuderink et al . found that neglecting the 
viscoelasticity of the tube wall could result in an 
underestimation of both phase velocity and damping. 
Hawley discovered that a viscoelastic wall model yielded 
results closer to normal physiological reality than an 
elastic wall model. 

In most computer models, arteries were treated as 
straight tubes. But the structure of the Circle of 
Willis is three-dimensional, and arteries are bent. 
Curvature together with pulsatile flow greatly affect 
wall shear stress and cause secondary flow which may be 
important in understanding the progression of 
neurovascular disease. Although curvature was not 
considered in prior intracranial cerebral circulation 
models, blood flow in extracranial curved arteries such 
as the aortic arch, carotid siphon and coronary arteries 
were analyzed by a number of researchers. Friedman and 
Ehrlich [Friedman and Ehrlich, Numerical simulation of 
aortic bifurcation flows: the effect of flow divider 
curvature , J". Biomechanics, 17:881-889 (1984)] digitized 
the radiography of human aortic bifurcation to obtain 
the contours of the artery and performed the computation 
of steady two-dimensional flow in the region. The 
results suggested that wall slope may be an important 
factor affecting the variability of shear stress along 
the medial wall. The variation in the curvature of the 
proximal iliac entries may affect the susceptibility of 
these arteries to vascular disease. Chang and Tarbell 
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[Chang and Tarbell, Numerical simulation of fully 
developed sinusoidal 

and pulsatile flow in the curved tubes , J. Fluid Mech. , 
161:175-198 (1985)] simulated pulsatile flow in the 
5 aortic arch and found that the secondary flow was nearly 
as large as the axial flow component, the secondary flow 
was complex with up to seven vortices and peak axial and 
highest r.m.s. wall shear stress were found at the 
inside wall. Also, the axial-flow direction was 
10 reversed at the inside wall. 

Ling and Atabek [Ling and Atabek, A numerical study 
^ of flow in curved tubes simulating coronary arteries , J. 

O Biomech. , 21:927-937 (1988)] simulated flow in the 

^ coronary artery using the same method of Chang and 

P 15 Tarbell. They found that flow was quasi-steady under 

resting conditions but markedly unsteady under excited 
Ul conditions. Only a single secondary flow vortex was 

Q found. Perktold et al . [Perktold et al , , Wall shear 

stress distribution in the human carotid siphon during 
p 20 pulsatile blood flow , J". Biomechanics, 21:663-671 

(1988)] used two different entrance velocity profiles to 
study the flow pattern in a slightly curved segment of 
the left main coronary artery. The different velocity 
profiles only resulted in significant differences in the 
25 immediate vicinity of the inlet, the difference does not 
persist far into the artery, and significant influence 
of the secondary flow on the wall shear stress 
distribution was found, 

Perktold et al . [Perktold et al . , Wall shear stress 
3 0 distribution in the human carotid siphon during 

pulsatile blood flow , J". Biomechanics r 21:663-671 
(1988)] analyzed flow in the carotid siphon and left 
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main coronary artery. They found that the maximum 
secondary flow velocities were on the order of three to 
four percent of the maximum axial velocity and the 
secondary flow has an important influence on the wall 
5 shear stress distribution. 

It is commonly believed that the blood behaves as a 
Newtonian fluid in the large blood vessel due to the 
high shear rate. However, some studies indicated that 
there were some non-Newtonian effects in the large 

10 vessel. It was also found that non-Newtonian effects 

exist in the low shear area near bends and bifurcations. 
Xu et al . [Xu et al . , Flow studies in canine artery 
bifurcation using a numerical simulation method , J , 
Biomech. Eng., 114:504-511 (1992)] used the Casson model 

15 to investigate the non-Newtonian effect in artery 
bifurcation. They found no great differences in 
velocity profiles. Lou and Yang [Lou and Yang, A 
computer simulation of the non-Newtonian blood flow at 
the aortic bifurcation , J". Biomechanics , 26:37-49 

20 (1993)] used a wea}c-form Casson model to study the non- 
Newtonian effect in aortic bifurcation. The results 
indicated that the non-Newtonian property of blood did 
not drastically change the flow pattern, but caused an 
appreciable increase in the shear stress and a slightly 

2 5 higher resistance to both flow separation and the phase 

shifts between flow layers. Dutta and Tarbell [Dutta 
and Tarbell, Influence of non-Newtonian behavior of 
blood flow in an elastic artery model , J". Biomech. Eng., 
118:111-119 (1996)] applied a simple power law model in 

3 0 an elastic artery and again found the viscoelasticity of 

the blood does not influence its flow behavior under 
physiological conditions in large arteries. 
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The present invention is a practical integration of 
several aspects of cerebrovascular and computer 
research. A program has been developed for computer 
aided neurovascular analysis and simulation that is 
5 useful for assessment and prediction of cerebral 

circulation changes. Such a program offers significant 
benefit, not only for stroke victims, but also for 
victims of trauma or congenital circulation disease. 
For instance in the case of progressive occlusion, the 
10 program offers the physician the ability to track 

circulation deterioration and plan for the eventual loss 
of the vessel by allowing the physician a practical 
^ method of evaluating potential treatment options . 

^ The program has up to four major components: (i) a 

O 15 vessel extraction system from Digital Subtraction 

is Angiography (DSA) ; (ii) a three-dimensional phase 

yi contrast Magnetic Resonance (MR) flow measurement system 

and three-dimensional pulsatility visualization; (iii) a 
vS computer simulation system for cerebral circulation; and 

- 

^ 2 0 (iv) a worldwide web-based interface and window to the 

world. 

The process involved in the present invention 
creates a virtual replica of the environs of the Circle 
of Willis or other vascular networks within the human 

2 5 body and computes the blood flow in this associated 

network. The modeling program also allows analysis and 
simulation of flow and pressure within the cerebrum 
and/or other organs under the condition of a selected 
blood flow perturbation (e.g., a cerebral aneurysm, 

3 0 stenosis, bypass, other cerebrovascular disease, trauma, 

etc.). The modeling program currently uses the finite- 
difference numerical method. Further under its various 
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embodiments, the process may be applied to any blood 
vessel leading to or coming out of the brain. Also, the 
same process can be applied to other vascular regions of 
the body . 

5 A vessel extraction system is used for determining 

exact vessel diameters for use in the model. The vessel 
extraction system may operate under any of a number of 
formats. For example under one format an Attentionally- 
Based Interaction Model may be used (e.g., "AIM: 
10 Attentionally Based Interaction Model for the 

Interpretation of Vascular Angiography" by F.K.H. Quek, 
C. Kirbas, and F. Charbel, IEEE Transactions on 
2 information Technology in Biomedicine, Vol. 03, No. 02, 

2C June 1999). Under another format (described below), a 

rj 15 semi-automatic process, involving the manual 

J . i 

2^ identification of the vessel wall, may also be used, 

U1 While either format may be appropriate, it has been 

found that the semi-automatic process described below 
J1 provides more accurate results . 

2 0 The vessel extraction system provides accurate and 
^^0 rapid measurement of vessel diameters and length in 

various predetermined regions of the Circle of Willis. 
The data with the smallest vessel size resolution (on 
the order of a few hundred microns) is currently derived 
25 from x-ray angiograms [platform: SGI Workstation, IRIS 
6.2/6.4 Varsity Package (Rapidapp, Open GL, Open 
Inventor, C++) ] . 

The use of digital x-ray angiography and 
computerized acquisition of the digitized picture can 

3 0 accelerate and enhance the data acquisition process over 

manual reading of vessel diameters from the x-ray 
angiogram. The two-dimensional nature of current x-ray 
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angiographic measurements limits, the ability to derive 
detailed three-dimensional data from the technique. 

Magnetic resonance imaging (MRI) has developed in 
such a manner that useful three-dimensional data 
5 pertaining to vessel size and flow are obtainable. The 
source of data is from magnetic resonance imaging which 
uses three-dimensional phase contrast angiographic 
methods to measure flow [platform: SGI Workstation, IRIS 
6.2/6.4 Varsity Package (Rapidapp, Open GL, Open 
10 Inventor, C++) ] . On its own, the vessel size resolution 
f the MRI data is on the order of a half a millimeter. 
The use of these data in the thre^dimensional 
phase contrast magnetic resonance (M^ flow measurement 
system of the invention enhances t^ accuracy of the 
15 cross-section and flow measurem^f^ts by three-dimensional 
localization and visualizatioj^of the vessels. Appendix 
A is a copy of the three-di^fensional phase contrast 
magnetic resonance f low nrfeasurement system and three- 
dimensional pulsatilit^ visualization system program. 
2 0 Turning now to image processing, the vessel 

extraction system will be discussed. For the 
interactive approach to be effective, it must exploit 
the facility of humans to operate or communicate within 
the bounds of certain rules which govern the approach. 
2 5 As will become apparent, the vessel extraction system is 
a selective interaction paradigm for the analysis of 
multimodal medical images (e.g., magnetic resonance 
angiography (MRA) , magnetic resonance imaging (MRI) , 
xenon computed tomography (XeCT) , x-ray angiograms 
30 (XRA) , etc . ) . 

FIG. 1 is a block diagram of a cerebral modeling 
system 10, generally in accordance with an illustrated 
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embodiment of the invention. The modeling system 10 
generally includes a data source, 12, a central 
processing unit (CPU) 14, a display 18 and a keyboard 
16. The data source 12 may one or more medical imaging 
5 systems (e.g., MRA, MRI, XeCT, XRA, etc.). 

In the system 10, the user selects the context from 
a menu (or a schematic diagram of the neurovascular 
system) and manipulates a cursor through the data using 
a 2D or 3D pointing device. For example, FIG. 2 shows 
10 the human Circle of Willis 100, including, 108 blood 

vessels. Selection of one of the vessels provides the 
system 10 with a context of what to look for . 
^0 FIG. 3 shows an interactive screen 110 that may be 

2r presented on the display 18 of the system 10. The 

Q 15 Circle of Willis 100 of FIG. 2 may be displayed as a 

graphical context representation in a box 112 of the 
in screen 110. An operator (not shown) of the system 10 

may select a vessel (e.g., 102) using the cursor 118. 

=.0 Once a context has been selected in the first box 

1 u 

1^ 20 112, the operator may show the CPU 16 where to look for 

'=0 the vessel in the image work area 114. The operator may 

identify the vessel 102 by placing the cursor 118 over 
the vessel and activate an ENTER button (click on the 
mouse) . 

2 5 However, as a first step, the user may also choose 

to calibrate the vessel extraction system using the 
screen 120 shown in FIG. 4. To calibrate the vessel 
extraction system, a known metal object (e.g., a penny) 
122 may be placed on the cheek of a patient during data 

3 0 acquisition. The vessel extraction system may be 

calibrated by measuring the metal object 122. 
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To calibrate the vessel extraction system, the user 
may place the cursor 118 on one side of the metal object 
122 and activate ENTER (click on the mouse) to identify 
a first side of the metal object 122. The user may then 
5 drag the cursor 118 across the metal object or simply 
click on the opposing side of the metal object 126. 

A distance measured by the vessel extraction system 
may appear in a first box 130. If the measured distance 
is accurate, the operator activates an ACCEPT button 
10 128. If the measured distance is not accurate, then the 
operator may adjust the pixel size to correct any error. 

Once the vessel extraction system has been 
calibrated, the operator may proceed to customize the 
p cerebral model to the precise characteristics of the 

p 15 patient. As a first step, the operator may select the 

^ screen 13 0 shown in FIG. 5 by selecting the appropriate 

Ln softkey along the lower edge of the control area 116. 

%, Shown on the right side in the image work area 114 may 

ffl be any niimber of different cranial slices. Shown in the 

ere : 

'^t 2 0 upper left in the graphical context box 112 is a portion 

^ of the cerebral vascular structure. 

Using a set of controls in the image control area 
116, the operator may view the entire vascular structure 
100, or any portion. The operator may zoom- in or out by 
25 selecting the appropriate softkey. 

Similarly, the operator may select any of a number 
of images for display in the image work area 114 . Shown 
along the left edge of the control area 116 is a set of 
softkeys for selecting the number of slices to be 
3 0 displayed in the work area 114. The CPU 14 

automatically adjusts a size of each frame to 
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accommodate the niimber of frames to be displayed in the 
image work area 114. 

To measure a vessel, the operator may select one of 
the frames within the image work area 114. To select an 
5 image, the operator may place the cursor 118 over an 

image and press ENTER. Alternatively, the operator may 
select the image by highlighting (i.e., clicking on) the 
image and selecting the softkey "MEASURE" in the left 
corner of the control area 116. 

10 Upon selecting an image, the screen 13 6 of FIG. 6 

may appear with the selected slice within the image 
display area 114. Softkeys in the process status area 
119 may be provided to allow the operator to zoom- in or 
out. By using the appropriate softkeys, the operator 

15 may position the appropriate vessel in a convenient 
location for measurement. 

To measure a vessel, the operator may first select 
the vessel to be measured in the context area 112 . 
Selecting the vessel may cause the selected vessel to 

2 0 appear in the color red. 

Upon selecting a vessel, the operator may then 
identify the boundaries of the corresponding vessel in 
the image work area 114. To identify the vessel 
boundaries, the operator may first place the cursor 118 

25 of a first edge 13 8 and activate the ENTER key. The 

operator may then place the cursor 118 or the opposing 
edge of the vessel and identify that edge to the CPU 14 
in a similar manner. 

The CPU 14 determines the distance between the two 

30 selected points 138, 140 and displays the measured value 
in a display box 142 . The operator may enter the value 
by selecting an accept box 144. Once entered, the 

28 



• 



measured dimension may be stored along with an 
identifier of the measured vessel. 

Once the operator has measured a corresponding 
vessel, he/she may return to the first screen 130 and 
5 select another frame and measure another vessel. The 
process of selecting frames and measuring vessels may 
continue until some or all of the vessels have been 
measured. 

Once the corresponding vessel has been identified 
10 and measured, a blood flow may be measured in that 

vessel. Blood flow may be determined using the three- 
^ dimensional location and orthogonal slice crossing that 

vessel. 

J The actual determination of blood flows is part of 

□ 15 a multi-step process referred to as cine phase contrast 

!Z magnetic resonance (PCMR) . Two-dimensional localization 

ill using cine PCMR requires positioning a line 

perpendicular to a centerline of a vessel of interest on 
0 a magnetic resonance angiography (MRA) image. Since it 

5i 20 is often difficult to position a line perpendicular to 

O the centerline, a three-dimensional process was 

de ve 1 oped . 

The 3-D localization method of slice selection was 
developed using a 1,5 T imager (manufactured by GE of 
25 Milwaukee, WI) . An axial MRA with zero interspacing 
covering the circle of Willis was gathered and 
transferred to a SGI workstation. A stack of 2D images 
can then be converted into 3-D volume data. The volume 
data can be described as follows: 

30 

V{x,y, z) =Ik(x,y) , 
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where k=0, 1, . . . N-1, z=k*T, I]^(x,y) is the kth MRA 

image, T is the slice thickness, and N is the number of 
slices. The x-axis is from right to left, y-axis is 
from anterior to posterior and z-axis is from inferior 
to superior . 

Given a volume of voxels V(x,y,z), an iso-surface 
in the form of a polygonal mesh can be generated using a 
Marching-Cube algorithm. For a given iso-surface value, 
each voxel can be classified into one of two types : one 
where the value at the voxel is greater than or equal to 
the iso-surface value. The other where the value is 
smaller than the iso-surface. 

According to the voxel type, a cube consisting of 
eight adjacent voxels may have 14 different 
configurations . The points of intersection between 
surface and the edges of the cubes can be generated by 
linear interpolation between values at the vertices. 
The triangular polygons may be used to tessellate these 
intersections . 

The iso-surface may be determined from a maximum 
intensity projection (MIP) using an appropriate 1.5 T 
imager (e.g., GE) as a data source 12. The procedure 
for determining such an iso-surface is as follows. 

The region encompassing the vessels of interest and 
the iso-surface value were obtained. The window level 
and the window width for changing the contrast of the 
MIP image were adjusted so that the vessels of interest 
were , enhanced and the background was suppressed. The 
iso-surface value was determined as follows: 

iso-surface value = window level- (window width/2) 
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Next, the iso-surface of cerebral vessels can be 
rendered by converting these triangular polygons into a 
scenegraph object using an appropriate software tool 
{e.g., "The Inventor Mentor: Programming Object-Oriented 
3D Graphics with Open Inventor, Release 2", Josie 
Wernecke, Addison-Wesley Publishing Co., 1994). The 
vessel of interest can easily be found through user 
interaction with the 3D scenegraph by rotating, 
translating, and zooming in and out on the 3D object. 
According to the orientation of the vessel, an image 
plane (either axial or sagittal or coronal cut) can be 
selected by the operator for localizing the vessel of 
interest. The plane 146 may be chosen such that it is 
visually closest to the plane that is perpendicular to 
the flow direction of the vessel (e.g., a sagittal cut 
for the middle cerebral artery shown in FIG. 9a) . 

FIG. 9a shows a sagittal plane for determining 
slice orientation for a middle cerebral artery. FIG. 9b 
shows the 3D vasculature obtained using the Marching- 
Cube Algorithm. 

Once the plane is selected, the 2D image associated 
with the plane can be calculated based upon the original 
volume of voxels and the vessel measured as described 
above. For an axial cut, the 2D image is the same as 
the original image. For a sagittal cut, the 2D image 
can be obtained through the expression: 

I(y,z)=V(X,y,z) , 

where X is the coordinate of axis x for the sagittal 
plane. For a coronal cut, the 2D image can be obtained 
through the expression: 



I(x, z)=V(x, Y,z) , 



where Y is the coordinate of axis y for the coronal 
plane . 

The cross-section of the vessel of interest in the 
2D image can then be specified by drawing a region of 
interest (ROD in the 2D image and the ROI can then be 
drawn in the 3D vasculature. The localization of the 
vessel of interest can then be achieved by correlating 
the ROI in the 2D image and the ROI in the 3D image. 

To obtain the border of the vessel in the 2D image, 
the ROI is refined by image segmentation using 
thresholding. The center of the refined ROI may then be 
calculated. By moving the plane along the vessel 
segment, a set of center points of the vessel may be 
identified. A line-fitting algorithm may then be used 
to fit these points to one line or more than one line. 
In the case of more than one line, the one with the 
maximum length and which line length exceeds a specified 
slice thickness may be chosen. The center of the line 
and an end point of the line may be saved and provided 
as an input to the imager 14 as the start and end 
location of the cine PCMR protocol in no particular 
localization mode. 

There may be many vessels including arteries and 
veins in a cine PCMR image. To identify those vessels 
in a cine PCMR image often becomes difficult when the 
slice orientation is obtained using an oblique cut. 
Identification of the vessel has been found to be much 
simpler if the slice orientation is obtained using the 
3D localization described above. The vessel of interest 



is usually proximate the center of the cine PCMR images 
when the 3D localization is used. The slice orientation 
can be drawn into the 3D vasculature when needed so the 
operator can see where the slice orientation and vessel 
of interest is actually located in the 3D vasculature. 
FIG. 9b shows the 3D vasculature obtained using the 
images shown in FIG. 7a. 

The binary masks corresponding to the ROIs drawn at 
each time point in the cardiac cycle for vessel diameter 
measurement may then be passed to the flow calculation 
program. FIG. 10 shows a screen 150 that may be 
presented to the operator . Shown in the image work area 
114 of the screen 150 is a plane passing through a 
vessel of interest. Shown in the graphical context area 
is a ROI 148, surrounding the vessel of interest. 

The operator may select a waveform softkey 152 to 
view bloodflow. In response, the screen 154 of FIG. 11 
may be presented to the operator. The FOI 148 continues 
to be displayed in the context area 112. However, in 
the image work area 114, a phase-contrast image of blood 
flow is now presented. 

Using the phase-contrast image, blood flow may be 
calculated based upon phase-shift information. The flow 
in a pixel may be calculated by multiplying the velocity 
of that pixel by the pixel area. The flow through the 
vessel may be determined by adding together the flows 
from all the pixels in the measured vessel. The cine 
PCMR may also record the ROI area, minimum, maximum and 
mean velocities in each ROI at each time point of the 
cardiac eye le . 

The flow waveforms and maximum velocity waveforms 
can then be plotted over the time course of the 
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heartbeat as shown in the work area 114 of screen 156 of 
FIG. 12. The vasculature and source of the bloodflow 
may also be retrieved by the operator as shown in the 
context area 112 of screen 158 of FIG. 13. 

Mean flow value provides a single quantity called 
volumetric flow rate (VFR) . The mean flow is calculated 
by summing up all the mean flow values (all of the VFRs) 
throughout the cardiac cycle and then dividing by the 
number of cardiac frames using the expression; 

SVFR„ 

VFR= — , 

N 



where VFR is the mean volumetric flow rate, n is the 
cardiac frame index and N is the total number of cardiac 
15 frames . 

As described above, vessel size and location may be 
obtained from MRA, MRI, XeCT or XRA data. The 
dimensional information derived from the data may be 
enhanced using MRI operating in the phase contrast mode. 
2 0 The use of phase contrast MRI allows volumetric blood 

flow to be determined by measuring blood velocity across 
a cross-section of each vessel. 

The enhancement of the accuracy and vessel size 
resolution arises from the three-dimensional 

2 5 reconstruction of the vessels that currently uses an 

interpolation scheme that is constrained by a piecewise 
smooth volumetric flow equation. The flow data are 
corroborated by transcranial Doppler measurements at 
predetermined vessel intervals. The flow velocity has 

3 0 also been found to vary as a result of vasoreactivity . 
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Benchmark values for cerebral autoregulation have been 
established and are continuously refined. 

The data analysis system supports both automatic 
and interactive extraction of the vessel cross-sections. 
5 A color coding scheme facilitates visualization and user 
interaction with the data to reduce the inter-user 
variability. 

For example, a user can view, move and freely 
rotate a three-dimensional picture of the 
10 cerebrovascular network 100. The user can also place a 
plane anywhere in the three-dimensional picture, then 
view the two-dimensional cross-section of the 
cerebrovascular network 100 at that location. 
^ The vessels are clearly discernible in the cross- 

O 15 section, with color coding to denote flow into and out 

^ of the planar cross-section. The pulsatility of blood 

In flow can be animated and visualized in the currently- 

developed user interface 150 in three dimensions. A 
vessel can be selected by the user for detailed and 
^5 20 graphical analysis of the flow through the vessel cross- 

■~n section that shows the pulsatile flow changes animated 

in time. 

The computer simulation system 10 for cerebral 
circulation, employs the outputs from both the vessel 

25 extraction system and the three-dimensional phase 
contrast MR flow measurement system to calibrate, 
customize and drive the cerebral circulation model 
[platform: PC (Pentium), Windows95/NT or DOS, Lahey 
Fortran 77]. The model is reconf igurable to account for 

3 0 person-to-person variability in the cerebrovascular 

network. The overall number of vascular segments in the 
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flow model can be increased or decreased as needed to 
form a customized model for each patient. 

The computer simulation system is flexible and 
accommodates empirical observations of measurements from 
the x-ray and magnetic resonance angiograms as well as 
from direct measurements of flow using flowmeters during 
surgery, and transcranial Doppler (TCD) at selected 
sites . 

The computer model is a one-dimensional, explicit, 
finite-difference algorithm based on a conservation of 
mass equation, a Navier-Stokes momentum equation, and an 
equation of state relating local pressure to local size 
of artery applied at each vessel segment [Kufahl & 
Clark, ASHE J. of Biomech. Eng. 107:112-122 (1985)]. 
Appendix Bi<s a copy of the current computer model 
prograij^i^ FIG. 29 is a flow chart of the modeling 
pro,0ram . ^ 

Since the arterial netwoi?Ks contain vessel loops 
{as well as many branchings^ the pressure and flow 
nodes are staggered thro^fghout the model.. Each vessel 
is divided into many >^gments ; the flow nodes are 
located at segment/ends , the pressure nodes at segment 
centers. Any mdJ.ti-vessel network configuration can be 
specified sa3?ely from the data file. Appendix C is a 
sample dara file with the computer parameters for a 
typicaic human model . 

The model is forced by one or more pressure or flow 
signatures at appropriate locations. Usually, a 
pressure- time signature at the root of the aorta, 
obtained from prototype measurements serves as the 
forcing function. 
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As shown (FIG. 29) , the model of the system 10 is 
first initialized 202 with initial pressures and flows 
at all nodes in all vessels. A cross-sectional area is 
calculated 202 for all points using the current 
5 pressure. The mass balance is determined 204 by 

calculating pressures at all points except the pressures 
at vessel junction centers. An inlet forcing function 
is invoked 206 to update pressure and flow sources. The 
internal junction boundary conditions may be evaluated 
10 2 08 by calculating flows at all junctions and the 
pressures at all junction centers. 

A momentum balance may then be determined 210 by 
calculating the flows at all points except the flows at 
junctions. A set of terminal boundary conditions may be 
.15 determined 212 by calculating the flows at the last 
nodes in all efferent vessels. 

Finally, the current time value is incremented 214 
and the incremented time is compared 216 with a modeling 
period. If the incremented time is less than the 
20 modeling period, the process 200 repeats. If not, the 
process terminates . 

A baseline vessel network 100 is used in the 
current model, including the Circle of Willis, 
ophthalmic arteries and other natural anastomoses . When 

2 5 surgical anastomoses are considered, the number of 

vessels can rise as needed. In addition, segments 
stenosed (perturbed) by any specified amount can be 
placed in any number of vessels. Aneurysms can also be 
simulated at various sites in the network. The results 

3 0 of any of a number of surgical procedures may be 

accurately predicted based upon the model results of the 
system 10. 
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A worldwide web-based interface and window to the 
world is available to provide blood flows and simulation 
results to physicians and surgeons of remotely located 
patients . 

5 Early models of blood flow in the brain were based 

upon an assumption that any particular artery feeds a 
particular part of the brain. The volume of brain fed 
can be calculated. First, the total mass of the brain 
can be determined. Then, the portion fed by a 
10 particular artery can be determined as a percentage of 
the determined total. Knowing the mass of brain fed by 
a particular artery allows the volume of blood necessary 
5 to feed that mass to be determined. Those early models 

J were general, using nominal values for the elasticity of 

Q 15 the blood vessels, the viscosity of the blood, and the 

^ vasculature arrangement of a normal patient's brain 

m circulatory system. 

The present invention is a refined model that is 
i,p capable of being adapted to specific patients. Once the 

20 volume flow of blood in certain arteries of a given 
patient has been determined, the model's terminal 
resistance pattern is determined for that patient. 

Deviations of the arterial structure of the blood 
supply of the patient's brain from the general model are 
25 identified from the angiograms. An x-ray angiogram (XR 
angiogram) of the patient's brain is used to determine 
the diameter of the blood vessels. Phase contrast 
Magnetic Resonance imaging angiography (MR angiography) 
is then used to determine an actual blood flow in 
30 critical arteries in the brain. Missing or additional 
arterial segments may be identified and used to adjust 
the model. A knowledge of the actual arterial structure 
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and actual blood flows can be used to customize the 
model to the actual patient. 

To validate the 3D localization of slice selection 
for cine PCMR imaging, a flow phantom was constructed 
shown in FIG. 14a. The cylinder is about the size of a 
human's head and filled with distilled water. Inside 
the cylinder there is a cube-like frame consisting of 
nine Plexiglas tubes of two different diameters. Those 
tubes are connected with each other sequentially. The 
tubes were circular in cross-section with inner 
diameters of 1/4 and 3/8 inches, i.e., 6.3 5 mm and 4.7 6 
mm. The fluid used in the experiment was 5 liters 
distilled water doped with gadolinium chloride-EDTA . 

To validate our flow measurement, both constant 
flow and pulsatile flow were evaluated. Axial 3D time 
of flight (TOF) images were acquired for the 3D 
localization. The images were extracted and transferred 
to a SGI workstation during scanning. . The 3D cube-like 
frame was then reconstructed using the marching cube 
algorithm described above. The coordinates required for 
an oblique cine PCMR were determined using the 3D 
localization. Cine PCMR images perpendicular to the 
tubes were acquired with those coordinates. A person's 
heart rate was used to trigger the sequence. All the 
images were extracted and transferred to a SGI OCTANE 
workstation for flow calculation. All the flows were 
evaluated using the software system CMIS implemented in 
C-f+, X/Motif, and Open Inventor on the SGI OCTANE 
workstation . 

To maintain a constant flow, the drained-out water 
was frequently poured back into the supply container. 
The flow rates were measured by timing the filling of a 
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graduated cylinder as the actual flow rates. The slice 
selections of the eight tubes were determined using the 
3D localization and validated with the spatial 
relationships of the tubes in the phantom. The cine 
5 PCMR were performed based on these slice selections. 

The PCMR flow rates were evaluated using the CMIS system 
and compared with the actual flow rate. Three different 
flow rates were achieved by adjusting the height of the 
supply container. The PCMR flow measurements were 

10 evaluated and also compared with the actual flow rates 
in the three different actual flow rates. 

Pulsatile flow was obtained using a valve system. 
An inline transonic flow probe was used for measuring 
the pulsatile flow rates for comparison. The PCMR flow 

15 measurements were evaluated and also compared with the 
manually-obtained flow measurements A misalignment of 
10°, 20°, and 30° for the central tube were calculated. 
The flows under the three different misalignments and 
the flow under no misalignment were evaluated and 

20 compared with the flow readings from the flow probe. 
The accuracy of the flow measurement under the three 
misalignments were obtained and compared with that under 
no misalignment. 

Two healthy volunteers were evaluated with cine 

2 5 PCMRA images using the 3D localization. A volunteer 

carotid study was performed to study the effect of the 
misalignment on the in vivo flow measurements. Three 
misalignments of 10°, 20°, and 30° for the left common 
carotid artery were calculated. The flows under three 

3 0 misalignments were compared with the flow measured under 

no misalignment. The second volunteer study was 
performed to validate the accuracy of the slice 
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selections for major arteries in the cerebral 
circulation determined by the proposed 3D localization. 

The slice selections for the major arteries in the 
cerebral circulation were calculated using the 3D 
5 localization. Cine PCMR were performed based on these 
slice orientations and the flows were evaluated. A 
peripheral gating was used to trigger the cine PCMR 
sequence. The following parameters were used in the 
cerebral blood flow protocol: TR, 33 ms; TE, minimum; 

10 flip angle, 30; number of excitations, 2; field of view, 
16 cm; image matrix, 256 x 128; number of phases, 32; 
slice thickness, 5 mm for carotids, velocity encoding 
(VENC) , 100 cm/sec for carotid, vertebral, and basilar 
artery, 70 cm/sec for other vessels. 

15 The 3D surface rendering of the phantom is shown in 

FIG. 14b. There are three lines in the center of the 
bounding box. The respective lines in FIG. 14b 
represent Right, Anterior and Superior directions 
respectively. The slice selections for eight tubes are 

2 0 shown in FIG. 15. The tubes were labeled from 1 to 8 

according to the slice selections shown in FIG. 15. 
Since the landmark of the head coil was set to the 
center of the cylinder, tube 1 is in the right and 
posterior position, tube 2 is in the superior and 
25 posterior position, and so on. The cine PCMR magnitude 
images for the eight tubes are shown in FIG. 16. The 
tubes of interest in FIG. 16 are all about in the 
centers of these cine PCMR images, circular and bright. 
The black rings shown in FIG. 16 are the tube walls. 

3 0 The large black areas and large gray areas shown in FIG. 

16 are outside the cylinder (i.e., air) and inside the 
cylinder (i.e., water), respectively. 
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The coordinates of the central points for all the 
eight slice selections were obtained from the 3D 
localization and given in Table 1. These coordinates 
are all consistent with the positions of tubes. Table 2 
gives the flow measurements calculated using the two 
methods described below from CMIS system and the 
accuracy of the flow measurements. The actual flow 
rates for all the tubes are the same since they are 
sequentially connected with each other. The actual flow 
rate is 33 0 ml/min. The difference between Method 1 and 
Method 2 is that Method 2 uses the flow offsets 
estimated from a cine PCMR imaging with no flow in the 
tubes while Method 1 does not use any flow offsets. The 
flow offsets of 16 ml/min and 8 ml/min for the tubes 
with diameters of 6.3 5 mm and 4.7 6 mm were calculated 
from the cine PCMR acquired using the same coordinates 
as tube 1 and tube 8, respectively, but with no flow in 
the tube . 

TABLE I show coordinates of the center points of 
the eight slice orientations. (R, A, and S stand for 
Right, Anterior, and Superior, respectively. The unit 
of the coordinates is mm. ) 



TABLE I 





1 


2 


3 


4 


5 


6 


7 


8 


R 


51.5 


8.1 


-49 .7 


-51.4 


-52 .7 


8 . 8 


47 . 5 


4 . 5 


A 


-49 . 5 


-45.2 


-51.7 


0 . 0 


51 . 5 


58 . 5 


54 . 0 


4 . 0 


S 


-5.3 


48.7 


2.7 


-62 .2 


-7.3 


34 .7 


-5.3 


-7.3 



TABLE II shows constant flow for eight tubes. (The 
actual flows for all the eight tubes are 33 0 Ml/min. 
Method 1 does not use any flow offsets while Method 2 
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uses flow offsets. The flow offsets for the tubes with 
diameters of 6.35 mm and 4.76 mm are 16 Ml/min and 8 
Ml/min, respectively. The units for flow are Ml/min, ) 



TABLE II 





1 


2 


3 


4 


5 


6 


7 


8 


Method 1 


347 .21 


320.03 


319.13 


337.89 


343 -77 


347.39 


318.82 


341.82 


Accuracy 


5.2% 


-3.0% 


-3.3% 


2.4% 


4.2% 


5.3% 


-3.4% 


3 .6% 


Method 2 


331.21 


336.03 


327.13 


329.89 


327 .77 


331.39 


326.82 


333 . 82 


Accuracy 


0.4% 


1.8% 


-0.9% 


0.0% 


-0.7% 


0.4% 


-1.0% 


1.2% 



TABLE III shows diameter measurements for eight 
tubes. (The unit for the diameter is mm) 



TABLE III 





1 


2 


3 


4 


5 


6 


7 


8 


Actual 


6.35 


6.35 


4. 67 


4. 67 


6.35 


6.35 


4 . 67 


4.67 


Measured 


6 .24 


6.24 


4.6 


4.78 


6.36 


6 .47 


4.6 


4.56 


Accuracy 


-1.7% 


-1.7% 


-1.5% 


2.3% 


0.2% 


1.9% 


-1.5% 


-2.3% 



FIG. 17 shows the comparison between the flow 
measurements for the tube located in the center of the 
phantom with diameters of 4.76 mm from cine PCMR and the 
actual flow rates in three different levels, i.e., 140, 
23 0, and 32 0 ml/min respectively. The average flow 
rates estimated from cine PCMR were 150.63 ml/min, 
237.70 ml/min, and 331.44 ml/min, respectively. 

FIG. 18 shows the comparison of pulsatile flow 
measurements for the tube located in the center of the 
phantom between PCMR and transonic flow probe. The 
average flow rates from both PCMR and the flow probe 
were 189 ml/min and 184 ml/min, respectively. 
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TABLE IV shows the effect of slice orientation on 
pulsatile flow measurements for the center tube. 
(Actual means the average flow rate or mean VFR obtained 
from transonic flow probe, MR means the mean VR obtained 
from PCMR measurement. The units are ml/min. ) 



TABLE IV 



Orientation 


Actual 


MR 


Accuracy 


0 


277 


284 


2.5% 


10 


282 


317 


12% 


20 


279 


341 


23% 


30 


283 


358 


26.5% 



FIG. 19 shows the 3D surface rendering of the first 
volunteer's carotid under different misalignments . 

TABLE V shows the effect of the misalignment on the 
flow measurements on the volunteer's left common carotid 
artery. (The units are ml/min.) 



TABLE V 



Orientation 


0 


10 


20 


30 


VFR 


545 


544 


514 


488 


Percentage 


0% 


0% 


-6% 


-10% 



TABLE VI shows the flow rates of major arteries in 
the cerebral circulation measured in the second 
volunteer study. (The units are ml/min.) 




TABLE VI 



LCCA 


LICA 


LVA 


LMl 


LM2 


Lpcom 


IiA2 


LPCA 


379 


226 


164 


176 


119 


21 


86 


81 


RCCA 


RICA 


RVA 


RMl 


RM2 


RAl 


RA2 


RPCA 


516 


381 


86 


197 


67 


196 


94 


88 



^0 



FIG. 2 0 shows some examples of the 3D localization 
obtained from the second volunteer study. The planes 
5 shown in FIGs . 12a, 12b, 12c and 12d give the slice 
selections for a left middle cerebral artery, a right 
anterior cerebral artery, a left posterior communicating 
artery, and a middle cerebral artery M3 branch, 
respectively . 

10 The eight planes shown in FIG. 15 represent the 

eight slice selections for the eight different tubes in 
the phantom. The center points of the intersections 
between the planes and the tubes are listed in TABLE I. 
The position of the first point (51.5, -49, -5.3) is at 

15 right and posterior area of the phantom. This verifies 
that the slice selection of the 3D localization is 
consistent with the location of the first tube. It can 
also be shown that the other seven slice selections are 
consistent with the positions of the seven tubes. In 

2 0 TABLE II, the flow measurements using Method 1 have less 
than ±5.3% error. Using Method 2, the accuracy is less 
than ±1.8% error. The diameter measurements were also 
given in TABLE III, The accuracy for the diameter 
measurements is less than ±2.3% error. 

2 5 To see the accuracy of flow measurement on the 

constant flow phantom under different flow rates, three 
different flow rates (140, 230, and 330 ml/min) were 
used. The flow measurements have an accuracy of less 
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than ±8.0% error without flow offset compensation. If 
the flow offset of 8 ml/min was considered, the flow 
measurements have an accuracy of less than ±2.0% error. 

In the pulsatile phantom study, the accuracy of the 
5 flow measurement in FIG. 18 has an accuracy of 2.7% 
error. The pulsatile flow waveform from PCMR is in 
excellent agreement with the flow waveform from the 
transonic flow probe. The effect of misalignment on the 
flow measurement of the pulsatile flow phantom can be 

10 seen in TABLE IV. The inaccuracy of the flow 

measurement is increased from 2.5% to 26.5% as the 
misalignment is increased from 0*^ to 30^. The 
misalignment effect in vivo is shown in TABLE V. Since 
there is no actual flow rate, the percentage changes of 

15 the flow measurements under misalignments with respect 
to the flow measurement under no misalignment were 
calculated. The percentage change of the flow 
measurement is increased from 0% to 10.0% as the 
misalignment is increased from 0° to 30^. 

20 The asymmetrical flow in the carotid arteries can 

be seen from TABLE VI in the second volunteer study. 
The explanation is that the volunteer had no left 
anterior cerebral artery and no right posterior 
communicating artery. This can be seen from FIG. 20. 

25 FIG. 20 shows not only that the slice selection can 
correctly be calculated, but also that PCMR flow 
quantification on current intracranial vessels can be 
achieved using the 3D localization. 

In summary, the 3D localization can be used in 

3 0 cerebral blood flow with PCMR imaging. The technique 
could prove to be of great value in both evaluating 
disease and increasing the effectiveness of clinical 
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applications of PCMR imaging. It also becomes possible 
for the hemodynamic parameters of some curved 
intracranial vessels in cerebral circulation to be 
evaluated through the non-invasive PCMR technique. 

An empirical study of user variability using the 
three-dimensional phase contrast MR angiographic flow 
measurement system was conducted. The study tested the 
repeatability of measurements by having ten users select 
the vessel cross-section from the interpolated color 
velocity images of the cerebrovascular flow. The flow 
computation performance of the system was verified 
against the data collected from a constant flow phantom, 
from ten patients with varying cerebrovascular disease, 
and from a computational fluid dynamic (CFD) simulation. 

The results of the empirical study show that the 
readings are stable across users, and that the flow 
measurements show a good degree of fidelity to the flow 
obtained through clinical measurement and CFD 
simulation. A 96-percent overall accuracy on 160 
measurements was observed for the phantom, with less 
than 5 percent mean error in the inter-user variability 
in extraction of the vessel cross-section among the ten 
users. For the ten patients, the flow computed was 
satisfactory in qualitative evaluation. The computed 
flow also correlated very well with those measured with 
a flowmeter in two patients r=0.985, P<0.001). The flow 
measurements in the Circle of Willis correlated well 
with the results of a CFD simulation of two patients (r- 
0.970, P<0.001) . 

The surgical treatment of giant aneurysms or head 
and neck neoplasms often requires permanent occlusion of 
the internal carotid artery. Tolerance to permanent 
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carotid artery sacrifice is traditionally assessed by 
temporarily occluding the carotid artery. This 
procedure is not without risk. The computational model 
of cerebral circulation of the present invention is a 
5 safer alternative to the balloon occlusion test (BOT) . 

The model used by the system 10 was used to create 
a virtual replica of the Circle of Willis and compute 
network blood flow using the finite-difference method. 
To evaluate the ability of the model to identify 
10 patients who tolerate permanent carotid occlusion, the 
difference in ipsilateral computed middle cerebral 
artery flow between patients passing and failing the 
balloon occlusion test was determined prospectively. 
:^ Each patient underwent four-vessel angiography and 

p 15 awake temporary occlusion of the internal carotid 

J£ artery. Failure of the BOT was defined as appearance of 

in any one of the following: a neurological deficit during 

the BOT, slowing of EEG wave patterns, and a drop in 
RSO2 of >10 percent. 
^ 2 0 Five of the 22 patients failed the BOT. The change 

C- in computed middle cerebral artery flow was 7 +/- 0.74 

cc/min for patients passing the BOT (p<0.001). 

The value of computational fluid dynamics (CFD) in 
predicting the outcome of BOT is striking. The results 
25 suggest a role for CFD in the evaluation of patients for 
permanent artery occlusion. 

The numerical model of the cerebral circulation was 
utilized along with the data from other supplemental 
diagnostic modalities to evaluate cerebrocirculatory 
3 0 collateral function during the balloon occlusion test 
(BOT) . 
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EEG, transcranial Doppler, cerebral oximetry, 
SPECT, were also done along with the computational 
modeling for patients undergoing temporary occlusion of 
the carotid to assess cerebral circulation when 
permanent occlusion is needed. 

Case example: A 49 year old female displayed 
diplopia and headache. Upon investigation, she was 
found to have a large right cavernous internal carotid 
artery (ICA) aneurysm. A computer analysis of her 
cerebral blood flow was made to assess her cerebral 
circulation, focusing mainly on the total middle 
cerebral artery flow. The numerical value of the total 
middle cerebral artery blood flow with and without 
occlusion of the ipsilateral ICA was predicted by the 
computer flow to be almost the same. 
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Vessel 


Flow computed before 
balloon occlusion 


Flow after 

balloon 
occlusion 


Left Middle 
cerebral artery 
and branches 


83 


81 


40 


40 


17 


16 


17 


16 


16 


17 


Total MCA Flow 


169 


172 ' 


Right Middle 
cerebral artery 
and branches 


80 


79 


38 


35 


18 


17 


17 


17 


17 


15 


Total RMCA Flow 


169 


172 



In view of the predicted toleration of the occlusion, 
the patient underwent balloon occlusion of the ICA on 
the right side, just distal to the aneurysms. The 
patient displayed good flow in both cerebral hemispheres 
with minimal changes in total computed middle cerebral 
artery flow. Concomitantly, there 'were no change in 
amplitude in EEG waves, recordings of cerebral oximetry 
using near infrared spectroscopy, transcranial Doppler 
(TCD) , or clinical deterioration. 




Temporary balloon occlusion is invasive and adds 
risk to the evaluation of patients who may already be at 
risk for infarction. The numerical model of the 
cerebral circulation is employed as an aid in the 
5 evaluation of patients for permanent internal carotid 
occlusion . 

Pulsatile pressure and flow in normal and diseased 
vessels (e.g., aneurysms, stenoses) can also be 
simulated. Modeling of network flow with anatomical 
10 variations, bypasses, or lepto-meningeal collateral 
vessels is also possible. 

In this particular case, this model was utilized to 
Q simulate changes in the middle cerebral artery blood 

"^-Z: flow before and after the balloon occlusion. If the 

C3 15 blood flow is decreased in a region of the Circle of 

T% Willis after occlusion of the carotid artery, the 

tfi simulation shows red at the region receiving decreased 

J' flow, and displays numbers corresponding to the 

O magnitude of the decrease. If the blood flow is 

iTi 20 increased in a region of the Circle of Willis after 

O occlusion, the simulation shows green at the region 

^ receiving increased flow, and denotes the magnitude. 

From the simulation output, the user can estimate 
whether the patient will tolerate the occlusion (or 
25 would pass the balloon occlusion test) . 

Under another embodiment of the invention, the 
modeling system may be extended to sectors of the 
cerebral vasculature. The generic model results, 
obtained using the individual patient's vessel sizes 
30 (obtained from digital subtraction angiography (DSA) ) 

and a generic pattern of peripheral resistance, indicate 
satisfactory accuracy for some, but not all, vessel 
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flows. In order to improve the simulation, patient- 
specific information is needed regarding the peripheral 
resistance distribution. To achieve this patient- 
specificity, the flows in certain crucial cerebral 
arteries must be obtained from phase contrast magnetic 
resonance angiography (PCMRA) . A patient-specific 
pattern of resistance can then be calculated. To this 
end, the original computer program may be modified so 
that all the efferents to a cerebral sector are melded 
into one terminal efferent vessel and a measured input 
flow to that sector is then used to determine its 
resistance. Seven such sectors have been identified in 
the cerebral macrovasculature . By combining the sector 
concept with PCMRA and DSA angiography, enhanced 
simulations of normal, diseased, and surgeon-modified 
cerebral configurations have been produced. 

Next we provide a generic model . The behavior of 
the cerebral arterial network with its pulsatile flow 
occurring in compliant arteries is mimicked using a 
computer algorithm for unsteady flow in a network of 
elastic vessels. The primary vessels of the cerebral 
macrovasculature are included in the network along with 
a skeleton of vessels, including the heart efferent, to 
represent the systemic vasculature. It was deemed 
necessary to include the aorta and subclavians in the 
simulation because these arteries are the staging 
vessels where important hemodynamic decisions are made 
(based on downstream blockage and resistance 
distributions) as to the delivery of blood into the 
carotids and vertebrals and the ultimate distribution of 
blood in the anterior and posterior sections of the 
circle of Willis. The circle acts as an invaluable 
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distribution center for cerebral blood but the 
distribution decisions are initiated in the aortic arch 
and subclavians. 

The numeric algorithm that is used to generate the 
pressures and flows at all locations within the vessel 
network and at all times during the heart pulse period 
is a one-dimensional, staggered, explicit finite- 
difference scheme, based on a conservation of mass 
relation to calculate the pressures, and a Navier-Stokes 
momentum equation to calculate the flows. As seen in 
the cross-hatched portion of the definition sketch below 
(FIG. 22), the advanced- time pressure can be calculated 
using the present-time pressure and the present-time 
spatial gradient of the flow according to the following 
difference equation : 

where Q is flow, P pressure, and A cross-sectional area. 

The compliance is given by the following logarithmic 
relation attributable to Raines: 



a(x, p) = Ao(x, pj 1 + — r In — 



(7) 



Aq (x, Pq ) Pq 

where p is a constant and subscript 0 refers to initial 
conditions. The underlying assumption of this relation 
is that the compliance is inversely proportional to the 

aA p 

pressure, i.e., — = — . 

dF P 

When the flow is in the positive x-direction, the 
advanced- time flow (see dashed portion of the definition 
. sketch above) can be calculated using the advanced- time 
spatial pressure gradient, the present- time spatial flow 




10 



15 



20 



25 



gradient, and the temporal flow gradient according to 
the following difference equation: 



qt-q: , k{q'J-K-JqiJ A'.Qr''A;.Qt 



Ax 



(8) 



where A and Q. are momentum and velocity shape factors 
respectively, and p, fx and ri are blood density, blood 
viscosity, and vessel inner radius, respectively. The 
last term of Equation 8, the resistance term, is 
evaluated using the slope of the velocity profile that 
is approximated by a sixth-order polynomial. 

Not all or even most of the arteries in the 
cerebral circulation can be included in the model. Some 
judicious choice must be made as to where to terminate 
the model- The terminal boundary condition imposed on 
all efferents of the chosen network is, therefore, most 
crucial to the flow-distributions predicted by the model 
and is the main thrust of this effort. A lumped 
parameter resistance-capacitance-resistance (RCR) 
combination was used in which the upstream resistance 
(Ri) was set equal to the characteristic impedance of 
the efferent. This relation when solved for the 
advanced-time terminal flow takes the following form 
(refer to RCR definition sketch) : 



pt+At 

^j-0.5 



Qr=- 



At R, 



(p;o.3-q;r.)^+5^ 



l + R, 



Cv^JL 

At R, 



(9) 



where Pj-0.5 and Qj are the pressure and flow at the inlet 
of the RCR, Cv is its volume capacitance, Pv is the 
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''venous pressure" at the outlet, Ri and R2 are the. 
resistances. 

The forcing function for the model is imposed at 
the root of the ascending aorta and is in the form of a 
5 normal pressure- time diagram there. A set of initial 

conditions, required to get the numeric process started, 
is given by setting all nodal pressures equal to the 
venous pressure and all nodal flows equal to zero. Ten 
to fifteen pulse periods of calculations are required to 
10 bring the pressures and flows into steady-state 

pulsatility. Computation time for a typical run on a 
400 MHz personal computer is approximately three 
g3 minutes . 

A sector model may be considered next. With the 
0 15 advent of PCMRA, the non- invasive measurement of 

H cerebral artery flow has become a reality. This advance 

p provided a means by which cerebral circulation 

simulations could be made more patient-specific. The 

0 crucial unknowns in such simulations are the terminal 

£1 20 resistances of the modeled vessels. If time in the MRI 

1 y 

p machine were not a factor, the flow in each individual 

efferent of the model could be measured and a 
corresponding resistance could be calculated. To limit 
machine time and the number of flows to be determined 

25 for the patient, a sector model of the cerebral 

circulation was devised in which all of the terminal 
vessels in any sector were melded into one terminal 
vessel. Seven such sectors are shown in FIG. 23 along 
with the primary and secondary flow inlets and/ or 

3 0 outlets and with the terminal efferent for each sector. 
To determine this one terminal resistance, it is 
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necessary to be able to define the flow into and/or out 
of the sector. 

This requirement takes different forms for 
different sectors. For the left middle cerebral sector 
5 (LBC) , there are three possible flow inlets (Ml, the 
main middle cerebral vessel, and two small secondary 
collaterals or anastomoses) and one terminal resistance 
efferent. This single efferent takes the place of four 
efferents used in the original generic model (FIG. 21) 
10 and has been formed by melding together these original 

efferents. The secondary anastomoses are not present iri 
all cases and, when present, it is difficult to measure 
O their flows. Therefore, the Ml flow represents the 

5 primary input to this sector and its mean flow value 

P 15 over a pulse period is acquired via PCMRA. The 

hj secondary flows are accounted for m the simulation by 

HJ means of their current calculated values as described 

5 below. The right middle cerebral sector (RMC) is 

^ handled in a similar manner. In the anterior cerebral 

pj 20 sector (AC), there are four possible flow inlets (the 

^ left and right anterior cerebral artery Al, and two 

\Q small secondary collaterals or anastomoses) and one 

terminal resistance efferent. In the left external 
carotid sector (LEG) , there are two possible flow inlets 
25 (the left external carotid artery and the left 

ophthalmic artery) and one terminal resistance efferent. 
The right external carotid sector (REG) is treated 
similarly. In the basilar sector (BS) , there are five 
possible flow inlets (the downstream segment of the 
3 0 basilar itself, the two posterior communicating 

arteries, and two small collaterals) and one terminal 
efferent (resulting from the melding of the two 
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posterior cerebral arteries and two superior cerebellar 
arteries). The basilar flow is the primary input to 
this sector and again it must be determined using PCMRA. 
The other four secondary vessels can be either inflow 
5 and/ or outflow vessels depending on the overall 

resistance balance in the system. The four secondary 
contributions are accounted for in the simulation by 
using their current calculated values. In the vertebral 
sector (VS) , there are three possible flow inlets (the 
10 left and right vertebral artery, and the downstream 
segment of the basilar) and one terminal resistance 
^ efferent. All three flows must be determined by PCMRA. 

start the computations for the sector model 
(FIG. 2\) , the assumption is made that the flow (QGOAL 
15 in FIG. 2X) in the terminal efferent for the sector is 

the same asNthe PCMRA measured flow in the primary input 
flow vessel, rgnoring, for this first set of 
calculations, tn^ other secondary inflows or outflows. 
A plausible set of\terminal resistances are needed to 

2 0 start the calculations and can be concocted either from 

a steady flow model of\;lie vessel system or from a 
previous similar simulatib^n. After each (pulse) period 
of calculation, the program ^ad justs the resistance for 
each sector terminal efferent ka^ multiplying its current 
25 resistance by the ratio of the current calculated 
terminal flow to the goal (PCMRA) fsLow for that 
terminal. This procedure would corre>sM:ly adjust any one 
terminal flow but, since all the other \erminal 
resistances are adjusted concurrently, errors are 

3 0 introduced and many adjustments are needed iti subsequent 

periods. This process continues for up to a u^r- 
specified number of pulse periods of calculation^^ 
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^s^od^max in FIG. 24) . A check is made at the end of 
each E5^iod to see if the measured sector flow goal has 
been achi^^d by the adjustment of the resistance of the 
terminal effeb^nt. If all seven sector flow goals are 
achieved within tnfe^period limit, the program proceeds 
to the next step. If Hae goals are not met after the 
maximum number of periods ,^^le program starts over with 
the current efferent resistancb^ and the original 
initial conditions of zero flow al:id venous pressure at 
all nodes in all vessels. 

The next step in the procedure is the adjustment of 
the flow goals (QGOALs in FIG. 24) so as to account for 
the secondary inflows and/or outflows. Since all 
vessels in the model are involved in the calculations 
discussed above, the secondary flows have been 
continually updated in this process. In this set of 
calculations, the flow in the terminal efferent is set 
equal to the primary PCMRA flow adjusted for the flows 
in the secondary inlets and/or outlets. For the left 
middle cerebral sector, this amounts to adjusting the Ml 
flow for the inflow or outflow contributions of the two 
secondary anastomoses. For the basilar sector, the 
contributions of the posterior communicating arteries 
and the two secondary anastomoses must be considered 
with proper signs. For the left external carotid 
sector, the ophthalmic flow must be combined with the 
left external carotid PCMRA flow in subsequent 
calculations. These adjustments are made only at the 
beginning of a period_max run, which can be cut short if 
all sector flow goals are met. When a user-specified 
number of adjustments have been made, the program ends. 
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At this stage of the simulation, a patient-specific 
set of resistances has been generated from patient DSA 
angiograms and measured PCMRA sector flows. This 
resistance pattern can now be assumed to be invariant in 
5 additional simulations of this patient. Modifications 
in system configuration can be made, say, by inserting a 
bypass vessel between the external and internal carotid 
arteries or by vessel occlusion in the treatment of an 
aneurysm. The neurosurgeon can now play ''what if?" 

10 games with the computer program to plan the best 
strategies for his/her treatment of the patient. 

Cerebral blood flow quantification may be 
considered next. Flow quantification with PCMRA 
provides a powerful tool to non-invasively analyze blood 

15 flow. It has been shown to do well in assessing both in 
vitro and in vivo flow velocities and volumetric flow 
rates. The inaccuracies of both velocity and VFR in 
vitro measurements in large vessels are about 5-10 
percent under continuous and pulsatile flow conditions. 

20 Blood flow measurements are influenced by many MR 

parameters: among these is the orientation of the 
vessel relative to the scan plane. The MR data 
acquisition is sensitive to any flow along the through- 
plane axis. If the through-plane axis is not parallel 

25 to the flux through the vessel, only a component of the 
true velocity is detected. This situation can arise 
when the imaging plane is not perpendicular to the 
vessel of interest. Wedding studied the effects of 
misalignment in vitro, 

30 A three-dimensional {3D) localization algorithm has 

been developed by Zhao et al . for the determination of 
slice selections that are perpendicular to the flow 
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directions for cerebral vessels during scanning and the 
identification of cerebral vessels from PCMRA images 
with oblique scans during post-processing. The 3D 
localization of slice selection was developed based on a 
1.5 Tesla MRI General Electric imager. A 3D time of 
flight MRA covering the circle of Willis is performed 
and the images are extracted and transferred to a SGI 
workstation. The 3D surface rendering of the 
vasculature is then reconstructed from these images 
using a Marching-Cube algorithm. The coordinates 
obtained from the 3D localization based on the 3D 
surface rendering specify the position of an oblique 
cine PCMRA scan. A cine PCMRA is then performed based 
on those coordinates. The following parameters were 
used: TR, 3 3 ms ; TE, minimum; flip angle, 30; number of 
excitations, 2; field of view, 16 cm; image matrix, 256 
X 12 8; number of phases, 32; slice thickness and 
velocity encoding (VENC) are 5 mm and 10 0 cm/ sec for 
carotids, VA, and BA, and 4 mm and 7 0 cm/ sec for MCA and 
ACA. The vessel area is semi-automatically determined 
through an interface on a SGI workstation. The velocity 
and volumetric flow rate (VFR) can then be calculated 
over the time course of the pulse period. 

This algorithm provides for the determination of 
cerebral blood flow in eleven major cerebral arteries 
from PCMRA (left and right CCA, left and right VA, left 
and right ICA, left and right MCA, left and right ACA, 
and BA) . The mean values of the VFRs from these eleven 
vessels provided the primary inlet/outlet flows used in 
the sector model. The MR scan time is about an hour per 
patient . 
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Following is a presentation of results and case 
studies. To exhibit the efficacy and abilities of the 
sector model of the cerebral circulation, three cases 
relating to the balloon occlusion test (BOT) used in the 
5 treatment of non-operable cerebral aneurysms will be 
presented and discussed. When an aneurysm is to be 
treated by occluding the vessel to which it is attached, 
a temporary occlusion of the vessel is made with a 
balloon to test the clinical outcome of the patient to 
10 permanent occlusion. Since the BOT is not without risk, 
the patient would benefit if a BOT computer simulation 
could be substituted. The three cases presented for 
discussion are intended to show how well the model 
1: predicts the BOT outcome. 

15 Case lo A 51-year-old white male presented with a 

W two-month's history of progressively worsening left 

pulsatile tinnitus. Angiogram and PCMRA studies showed 
^ that both the left and right posterior communicating 

S arteries and the Al segment of the left anterior 

5jf 20 cerebral artery were small. The left internal carotid 

artery contained a large aneurysm and a stenosis of 6 0 
percent. These data were used as Baseline in the 
computer simulation. The patient subsequently passed a 
BOT. A permanent occlusion of the LICA was made with 
25 balloons trapping the diseased distal cervical segment 
of the LICA. 

Case 2o A 56-year-old black female presented with 
a subarachnoid hemorrhage. Angiogram and PCMRA studies 
showed multiple aneurysms including a giant aneurysm at 
30 the left superior hypophyseal artery as well as a 40 
percent stenosis in the left internal carotid artery. 
The patient failed the BOT under hypotensive challenge. 
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becoming hemiparetic and aphasic during the test. SPECT 
also showed a massive difference between the two 
hemispheres with a significant decreased relative blood 
flow in the left hemisphere involving the 
5 frontoparietal, occipital and temporal lobes. A PCMRA 
CBF study was performed before the BOT. These data were 
used in the Baseline simulation. 

Case 3o A 71-year-old white female presented with 
persistent headache. Angiogram and PCMRA studies showed 

10 a large right carotid-ophthalmic aneurysm as well as 
occlusions of both left and right posterior 
communicating arteries, occlusion of the anterior 
communicating artery, and a 40 percent stenosis in the 
left anterior cerebral artery. These data were 

15 considered as the Baseline in the computer simulation. 

A PCMRA CBF study was performed before the BOT. During 
the BOT in the right internal carotid artery, an 
immediate failure ensued when left side hemiparesis 
developed . 

2 0 A typical angiogram, 3D localizer image, and flow 

waveform are presented in FIG. 25. The PCMRA and 
simulation results are shown in the form of bar graphs 
in FIG. 26 for Case 1 and in FIG. 27 for Cases 2 and 3. 
A typical set of simulation results is shown in both 

2 5 numeric and color forms for Case 2 in FIG. 28, These 

results will be discussed in detail in the next section. 

The purpose of this discussion is to demonstrate 
the efficacy of a novel simulation technique, the sector 
model, in the patient-specific prediction of clinical 

3 0 outcomes accompanying surgical intervention. The 

specific intervention considered here is the treatment 
of inoperable aneurysms by permanent occlusion of a 
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carotid artery. The patient's tolerance of such an 
occlusion is determined by temporarily occluding the 
artery, the balloon occlusion test (BOT) . . Since this 
test is not without risk, it would be beneficial for the 
5 patient to substitute a valid computer simulation for 
it. 

The PCMRA flow measurements along with the sector 
and generic model results are displayed in FIGs . 2 6 and 
27. For case 1, the Baseline results in FIG. 26a show 
10 excellent agreement between the PCMRA flow and the 

sector model results for most vessels. This result was 
expected since the PCMRA flow measurements were used in 
^2 the RCR adjustments to satisfy the iterative process of 

=p FIG. 24. Once the RCR pattern had been determined for 

2 15 this patient, it was used independently to predict the 

IJ post-BOT results shown in FIG. 26b. Here again, the 

agreement between the PCMRA and the sector model results 
B is excellent (with the possible exception of the LVA) . 

^ The generic model compared well with the PCMRA 

fy 2 0 measurements for some vessels, but not as consistently 

^ as the sector model. For cerebral perfusion pressure 

^0 (CPP) results in FIGs. 26c and 26d, both sector and 

generic models give comparable predictions for Baseline 
and BOT with no dramatic pressure drops during the test. 
25. Since the patient passed the test, these results are 
also compatible with the clinical outcome. 

Similar results were obtained for Cases 2 and 3 for 
Baseline CBF as shown in FIGs. 27a and 27d. Since both 
patients failed the BOT test, there are no post-BOT 
3 0 PCMRA data to compare with the simulation results. 

Therefore, only calculated Baseline and BOT CBF values 
are shown in FIGs. 27b and 27c for Cases 2 and 3, 
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respectively. These cases have drastically different 
flow patterns in the carotid vessels (RCCA, . LCCA, RICA, 
LICA) because of the balloon closure of the LICA and 
RICA in Cases 2 and 3, respectively. For the BOT 
5 simulation, the reduction in flow values for the LMCA in 
Case 2 and the RMCA for Case 3 in comparison with the 
Baseline values is compatible with the clinical outcome 
of BOT failure. A similar result for CPP values for 
these arteries is shown in FIGs . 27c and 27f. The 

10 reduction in flow and pressure for Case 2 was not as 

pronounced as that for Case 3 because Case 2 failed the 
BOT only under hypotensive challenge. In Case 3, the 
reduction in flow is nearly 100 percent and the pressure 
is reduced to venous pressure levels in the RMCA, 

15 For Case 3 , when the seven-sector model was used to 

simulate the BOT, an anomaly was discovered in the AC 
sector. All possible inflow paths to the RMC and -AC. 
sectors were seemingly blocked by the occlusions 
indicated in Section 5. However, the sector model found 

2 0 the one unblocked path. That path started by following 
the LAC up to the junction of the terminal efferent. 
Since only part of the available LAC flow was needed to 
satisfy that efferent, the rest passed down the RAC and 
went on to satisfy the RMC sector. To correct this 

25 implausible solution and to prevent the anomaly, an 
eight-sector model was used with two separate AC 
sectors. However, the resulting RAC and RMC sectors 
had no connection to the forcing function in the 
ascending aorta except for minor contributions from the 

30 small ophthalmic and two secondary anastomoses. the 
results of the new simulation showed essentially zero 
flow and venous pressure levels in these sectors. The 
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immediate left side hemiparesis in this patient is 
subsequently explained by these simulation results. 

An alternative method of presenting the simulation 
results is shown in FIG. 2 8 where the arterial network 
5 is depicted with the various arteries colorized (gray 
scales in printed version) in accord with the amount of 
flow in them. The actual flow values are also provided 
along side of the artery. The color coding agrees with 
that used in other modalities such as SPECT in order to 
10 provide the surgeon with an overall visual understanding 
of the flow changes resulting from a surgical 
intervention, in this case a carotid occlusion during 
BOT. 

In conclusion, patient-specificity has been 
15 achieved in the simulation and analysis of the three BOT 
patients. These results demonstrate a significant 
advance in computer modeling of the cerebral circulation 
through the introduction of the sector model. While the 
present study only examined three cases, the results 
2 0 clearly demonstrate the working hypothesis of this 
sector model. Many improvements are envisioned for 
future models: the forcing function can be made more 
patient-specific; the systemic portions of the model can 
be made into sectors and governed by PCMRA measurements; 

2 5 and the autoregulatory aspects of the cerebral 

arterioles can be added to the simulation. In addition, 
the model can be generalized to other regions of the 
cardiovascular system like the coronary circulation or 
the circulation in the extremities . A model that can 

3 0 accurately predict cerebral blood flow will provide an 

invaluable tool for physicians in managing complex 
cerebral problems . 



65 



The evolution of models of the cerebral circulation 
has paralleled the progress in the fields of fluid 
dynamics and computer science. With each new 
possibility of simulation, a closer approach to 
"physiological" modeling was achieved. The ideal model 
would be patient specific, highly predictive, and 
reflective of actual conditions both at the macro and 
micro circulation level. In addition, such an ideal 
model would be easily re-configurable to conform with 
the "real time" demands of clinical medicine. 

Although computer models of various degrees of 
sophistication have been proposed to simulate the 
cerebral circulation, they have been used mainly as 
theoretical and research tools. Turning them into tools 
for clinical applications has the potentially powerful 
benefit of predicting the result of neurovascular 
reconstructive procedures. However, none of the 
published models has ever been validated in vivo using 
quantitative h\iman data, partly because of the lack of 
available tools for measuring such data. 

The method described above offers significant 
benefit, not only for stroke victims, but also for 
victims of trauma or congenital circulatory problems. 
The model has also been developed for other parts of the 
body such as the heart, extremities, etc. Also it may 
have application to forensic uses. In any application, 
it is believed that the model offers significant 
benefits in helping to understand the effects of 
potential treatment options available to the physician. 

A specific embodiment of a method and apparatus for 
modeling cerebral circulation according to the present 
invention has been described for the purpose of 



illustrating the manner in which the invention is made 
and used. It should be understood that the 
implementation of other variations and modifications of 
the invention and its various aspects will be apparent 
to one skilled in the art, and that the invention is not 
limited by the specific embodiments described. 
Therefore, it is contemplated to cover the present 
invention any and all modifications, variations, or 
equivalents that fall within the true spirit and scope 
of the basic underlying principles disclosed and claimed 
herein. 



